% ---------------------------------------------------------------------------
%   2-D, MULTIGROUP DIFFUSION-BASED, NONLINEAR RESPONSE MATRIX SOLVER - TEST!
% ---------------------------------------------------------------------------
%  j. roberts, 01/23/2010
%   This MAIN file currently only drives a power iteration like scheme.  As
%   soon as possible, the nonlinear solver will be used.  This file depends
%   on:
%       redblack2D.m      - iterative solver
%       resp2dexample.m   - response function generatory (i.e. a diff. code)
%       twoDcoefMGresp.m  - coefficients for diff solve
%       cresp.m           - current responses
% ---------------------------------------------------------------------------
%profile on

clear;% clc;
format long
global order numg numel

%----------------------- GET PROBLEM DATA -------------------------
disp('----RB----')
tic
for or = 6:6% % .02457/0.02481
order = or; numg = 2; numel = 1;
j = ones(1,numel*4*(order+1)*numg); % a uniform guess of coefficients
% j = zeros(1,numel*4*(order+1)*numg); 
% for gg = 1:numg
%     j( gg:(order+1)*numg:end)=1;
% end
j = j / sqrt(j*j');
if numel == 1
    kk= 0.06892; % reference keff
    kk2=1.099009048640218;
else
    kk=1.096296747503673;%1.139551748660189;
    kk2=1.096561615874126;
end
x = [j .07]';
tol = [1.e-8 1.e-8]; % tolerance for norms
[x1,normk,it,nrm] = redblack2D(x,12);
disp([num2str(or),'  ', num2str( x(end) ), '   ', ...
        num2str(100*(kk-x(end))/x(end)), '    ', ...
        num2str(100*(kk2-x(end))/x(end)) ])

end
disp(['rb keff = ',num2str(x(end))])
t1=toc;
disp(['rb time= ',num2str(t1)])
return
%---------------------- NEWTON SOLVER -----------------------------
% uses the standard newton method w/ a fd-approximated jacobian
disp('----NEWTON----')
gm = 0;
tic
x3          = [x' 1.0]'; % add in current eigenvalue
xx3(:,1)    = x3;
z           = respfct2Dexample(x3); % initial residual
itmx        = 20; 
it          = 1;
zz(1)       = norm(z);
nrmcn(1)    = 1;
while ( zz(it) > tol(1)*zz(1)+tol(2) && it < itmx )
    fp          = jacob('respfct2Dexample',x3,z);
    s           = -fp\z;
    x3          = x3 + s;
    it          = it+1; 
    xx3(:,it)   = x3; % Keep all x's to estimate rho
    z           = respfct2Dexample(x3);
    zz(it)      = norm(z);
    x3(end-1)
    disp(zz(it))
end
t2=toc;
disp(['nm time= ',num2str(t2)])

profile viewer
return
